Calibration of microphone arrays with an uncalibrated source

ABSTRACT

Microphone array calibration that does not require a calibrated source or calibrated reference microphone is provided. We provide a statistical (Bayesian) algorithm that (under condition of reasonable environment noise during calibration) can determine gain and phase differences of a whole array at once, even when the gain and/or phase of the source is unknown. More specifically, a Bayesian regression with complex log-normal prior and complex normal likelihood is employed. The inherent phase-wrapping ambiguity in this regression is resolved by exploiting the similarity of likelihood between a lattice point and its Euclidean Voronoi region.

CROSS REFERENCE TO RELATED APPLICATIONS

This application is a continuation of Ser. No. 16/244,918 filed on Jan.1, 20119 and hereby incorporated by reference in its entirety.

Application Ser. No. 16/244,918 claims the benefit of U.S. provisionalpatent application 62/616,884, filed on Jan. 12, 2018, and herebyincorporated by reference in its entirety.

FIELD OF THE INVENTION

This invention relates to improved calibration of microphone arrays,e.g. by providing calibration without any need for a calibrated sourceor calibrated reference sensor.

BACKGROUND

Over the past decades, noise pollution has become an increasing problem.To reduce noise emissions, it is imperative that engineers have easyaccess to sound measurement equipment so that acoustic product designcan be improved. To this end, equipment such as a 1024 channel MEMSmicrophone array for near-field acoustical holography and far-fieldbeamforming applications is commercially available.

To get good near-field acoustical holography and far-field beamformingresults, it is important to compensate for differences between gain andphase deviations of the individual microphones. The goal of calibrationis to reduce the uncertainty of these differences. Conventionally thiscalibration is done using a set of measurements with known source.

SUMMARY

However, we want customers to be able to calibrate their microphonearray ‘at home’/in their office, i.e. without access to expensiveacoustic laboratory equipment such as a calibrated acoustic source. Soour source is fully or partially unknown. This is an aspect not coveredby existing microphone array calibration methods, to our knowledge. Moregenerally, there is no need for any calibrated reference in our approachfor providing calibration of arrays of acoustic microphones.

We provide a statistical (Bayesian) algorithm that (under condition ofreasonable environment noise during calibration) can determine gain andphase differences of a whole array at once, even when the gain and/orphase of the source is unknown. This guarantees a compensated array willalways outperform a (non-compensated) factory array. As part of thismethod, the phase-wrapping ambiguity is dealt with via a novel approach.

This is possible by taking into account the uncertainties in acousticsource, waveguide and the microphones (i.e. specification from themanufacturer or results from a prior calibration). Because there is aplurality of microphones, the (previously unknown) gain and phase of thesingle source can be estimated. These estimated properties can then inturn be used to calibrate the microphones. More technically, a presentlypreferred algorithm implements a Bayesian regression with complexlog-normal prior and complex log-normal likelihood. The inherentphase-wrapping ambiguity in this regression is resolved by exploitingthe similarity of likelihood between a lattice point and its EuclideanVoronoi region.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 schematically shows an acoustic microphone array calibrationsetup.

FIG. 2 show steps of a method for microphone array calibration accordingto an embodiment of the invention.

FIG. 3 shows steps of a phase unwrapping method suitable for use inembodiments of the invention.

FIGS. 4A-D are sketches corresponding to the steps of the method of FIG.3.

FIG. 5A shows a first exemplary acoustic waveguide configuration.

FIG. 5B shows a second exemplary acoustic waveguide configuration.

DETAILED DESCRIPTION

A) General Principles

Microphone array calibration is a topic of general interest: awell-calibrated array produces less disturbance, and reconstructioncould be improved by accounting for any remaining calibration errors.Array calibration has been studied before in the contexts of radar andbeamforming, but to the authors' knowledge no calibration method isavailable that:

1) can correct both gain and phase deviations of the sensors;

2) can guarantee a calibrated array will always outperform a(non-calibrated) factory array; and

3) does not require a precise reference signal (e.g. from an acousticsource or a reference sensor).

In this work we provide a novel calibration method with these propertiesin a Bayesian framework.

To begin, it will be helpful to consider a microphone array calibrationsetup as shown on FIG. 1. In this example, a single acoustic source 102provides a source signal s to a transmission medium 104 having transferfunction w from the source to each element of the microphone array 106.It is convenient to combine the source signal s and the transferfunction w into a parameter t which is the input to array 106. Hereboldface quantities denote N dimensional vectors where N is the numberof elements in the microphone array. The measured output of themicrophone array is the observation o. The quantity x accounts for thegain and phase variation of the elements of microphone array 106.Microphone array calibration amounts to providing an estimate for thegains and phases x based on the observation o and on the transferfunction w. For Bayesian estimation, an estimate for x prior tocalibration is also used.

FIG. 2 shows steps of a method for microphone array calibrationaccording to an embodiment of the invention. Step 202 is providing anacoustic source. Importantly, there is no requirement that this sourcebe calibrated. Step 204 is providing an estimate of the transferfunction (w) from the source to each array element of the N-elementmicrophone array. As seen below, it will suffice for this estimate to bea probabilistic estimate. Step 206 is providing measurements (o) fromthe array elements when the source is operating. Step 208 is performingBayesian inference of gains and phases (x) of the microphone array basedon the measurements and on the estimate of the transfer function and onthe estimate of prior x.

As used herein, Bayesian inference is defined as a method of statisticalinference in which Bayes' rule (i.e., P(A|B)=P(B|A)P(A)/P(B) where A andB are events) is used to find a (posterior to calibration) estimate ofthe gains and phases (x), based on beliefs from one or more (possiblyuncertain) measurements, preferably constrained with an informative(i.e. sufficiently known to reduce uncertainty in the measurement(s))belief about the gains and phases (prior to calibration). The beliefsare quantified as probability distributions. Thus, a quantitative methodis obtained to reduce uncertainty in the measurements (e.g., a fully orpartially unknown source, or deviations in the transmission medium, ornoise) with informative prior beliefs, thus improving the quality ofcalibration over the situation where only the measurements (and not theprior beliefs) would have been used. This is especially relevant forarray calibration, because the plurality of microphones amplifies theuncertainty-reducing effects of (possibly mild) prior beliefs about eachsingle microphone and hence can significantly improve quality of thecalibration. Note that if the measurements are not constrained withprior beliefs, or the prior beliefs are not informative, at least somecertainty about the measurements (e.g. source and waveguide) must beprovided to obtain usable calibration results. However, in the contextof array calibration, sufficient prior beliefs are typically readilyavailable in practice, e.g. from specifications of the microphonemanufacturer. The following description provides an illustrativedetailed example of a presently preferred approach for such Bayesianinference as applied to microphone array calibration.

In this example, the result for gain is a normal distribution and theresult for phase is an infinite weighted sum of normal distributionswith weights γ(k). Here k is an N-dimensional vector of integers. Step210 is phase unwrapping (i.e., truncating the infinite sum to a finitesum) by sampling a probability distribution function (PDF) of γ(k) andselecting the K best k values from that sample set. This approach forphase unwrapping is called shotgun unwrapping.

FIG. 3 shows steps of a phase unwrapping method suitable for use inembodiments of the invention. This method shows the sub-parts of step210 on FIG. 2 in greater detail. Here step 302 is sampling from acontinuous PDF of γ(k) to provide an initial k-set

₁. This can readily be done, since as seen below γ(k) is normallydistributed with known mean and covariance matrix. Thus standardstatistical methods can efficiently provide this sampling. Step 304 isto round the elements of

₁ to the nearest integers (more precisely, to the nearest lattice pointsin N-dimensional space) and eliminate any resulting duplicates toprovide a discretized k-set

₂. Step 306 is to evaluate the distance (using an appropriate metric asshown in detail below) between each element of

₂ and the mean of the PDF of γ(k). Smaller distances correspond to morelikely weights. Thus step 308 is to select the K elements of

₂ having the shortest distances as the K best k values. Here K is apredetermined integer that can be selected based on practicalconsiderations. K is limited by computational resources (i.e. larger Kresults in more memory requirements and longer runtime). Also, there isno need to make K larger than the amount of elements remaining in

₂. In some situations encountered in practice, the amount of remainingpoints in

₂ is already quite small, which makes selection of K trivial.

Selecting the K smallest distances from

₂ can be expedited according to known principles, such as removingelements of

₂ having distances greater than a predetermined threshold prior toselecting the K best weights. This is helpful since sorting

₂ is typically done as part of selecting the K smallest distances from

₂, and the sorting will take less time if the k values having distancesthat are way too high are removed first. In preferred embodiments, thedistances are computed as follows. Let the probability distribution ofγ(k) have a mean μ_(u) and a covariance matrix Σ_(u). Then the distancesM(k) are given by M(k)=√{square root over ((k−μ_(u))^(H)Σ_(u)⁻¹(k−μ_(u)))}.

FIGS. 4A-D are sketches corresponding to the steps of the method of FIG.3. Here FIG. 4A corresponds to step 302, where the sample points areblack disks and the mean of γ(k) is an open triangle. FIG. 4B shows theresult of step 304. Note that sample points are now only present atlattice points (i.e., intersections of the grid lines). FIG. 4C showsthe result of step 306. Every sample point has its distance to the mean(dashed line) calculated. FIG. 4D shows the result of step 308. Only theK sample points closest to the mean are retained (here K=5 to provide asimple example). Note that it is possible (as seen on FIG. 4D) for anearby lattice point to be missed by this method, but probabilisticallyit provides good results.

Although the example of FIGS. 4A-D appears to be trivial, since thereare much simpler and more direct methods to find lattice points close tothe triangle on these figures, that is an erroneous impression caused bythis example having N=2 for ease of illustration. As indicated below,finding these closest lattice points is NP-hard in the number ofdimensions N. Since relevant acoustic microphone arrays can havehundreds of elements or more, it is important to have a phase unwrappingapproach that scales well to high dimensionality.

As indicated above, an important advantage of this approach is that thesource doesn't need to be calibrated. The amplitude and phase of theacoustic source can be assumed to be drawn from a predetermined sourceprobability distribution. Such a source probability distribution can besignificantly ill-defined (e.g., unknown phase, amplitude range=dynamicrange of the source) without significantly impairing the calibration.This allows one to calibrate an acoustic microphone array with readilyavailable sources, such as the speaker of a smart phone, or moregenerally from any mobile electronic device having a speaker. Examplesinclude: An acoustic calibrator (or pistonphone) with known gain andfrequency but unknown phase. A smart phone or computer loudspeaker withunknown gain, unknown phase, and possibly unknown/unstable frequency.

Practice of the invention does not depend critically on the details oftransmission medium 104 on FIG. 1. As indicated below, this transmissionmedium can even be free space (preferably in an anechoic room). In somecases it is preferred to provide greater control over the transferfunction w by providing an acoustic waveguide network configured tocouple the acoustic source to the array of acoustic microphones. Theacoustic waveguide network can include a source port (e.g., 508 on FIG.5A) corresponding to the acoustic source, and array ports (e.g., 510 onFIG. 5A), each array port corresponding to a corresponding one of theelements of the array of acoustic microphones.

FIG. 5A shows a first exemplary acoustic waveguide configuration. Inthis example, acoustic source 102 is coupled to 1-D microphone array 506via acoustic waveguide network 502. Acoustic waveguide network 502 canbe implemented as a tree-like network of tubes 504. This amounts to a1×5 acoustic splitter. Here 508 is the source port and 510 are the arrayports of the acoustic waveguide network. FIG. 5B shows a secondexemplary acoustic waveguide configuration. Here a 2-D array ofmicrophones 540 is coupled to acoustic source 102 via an acousticwaveguide network including 1×5 acoustic splitters 512, 522, 524, 526,528, 530. For simplicity of illustration, coupling is shown on FIG. 5Bwith single lines, and lines that cross an element of array 540 aren'tcoupled to that element.

Acoustic waveguide networks can be fabricated via rapid prototyping(e.g., 3D printing). An advantage of 3D printing is that customers withaccess to a 3D printer could download their own waveguide design andfabricate it on-site without a manufacturer needing to stock and ship itto them.

Optionally, a calibrated reference microphone (e.g. IEC 61672 soundpressure level meter) can be used to improve the calibration. When thismicrophone has tighter manufacturing/calibration tolerances than themicrophones in the array, the calibration results can be furtherimproved. Because this calibrated microphone is allowed in formal soundpressure level measurements (e.g. as evidence in a lawsuit),incorporating such an ‘official’ microphone in the calibration allowsusers to make measurements from our microphone array moreaccepted/traceable for formal purposes. Such a calibrated referencemicrophone can be expensive and hence not always available ‘at home’.The main advantage of the present approach is that, unlike many existingmethods, the reference microphone is not necessary for the calibrationprocedure to succeed.

B) Mathematical Development

Notation and variables are defined in the tables at the end of thissection.

B1) Definitions

B1a) Complex Normal Distribution

If a is a complex normal random variable, then:a˜

(a; μ _(a), Σ_(a))=π^(−N) det(Σ_(a))⁻¹ exp(−∥a−μ _(a)∥_(Σ) _(a) ²),  (1)where:∥a−μ _(a)∥_(Σ) _(a) ²

(a−μ _(a))^(H)Σ_(a) ⁻¹(a−μ _(a))  (2)is the squared Mahalanobis distance.B1b) Complex Log-Normal DistributionIf:a˜

(a; μ _(a), Σ_(a)) and b˜

(b; μ _(b), Σ_(b)),  (3)then c is a complex log-normal variable:c=exp(a+ib)˜

(c; μ _(a) +iμ _(b), Σ_(a) +iΣ _(b)).  (4)B2) Algorithm DerivationFor measurements over a frequency range, the observation can betransformed to the frequency domain and the steps below can be performedfor each frequency bin of interest.Step aModel the observation as a multivariate complex normal distribution:o|s, w, x˜

(o; s·w⊙x, σ ² I)  (5)The location parameters (i.e. means) are set to the expected sound levelof the source times the transfer functions (from source to each sensor)of the waveguide. The scale parameters (i.e. covariance matrix) accountfor measurement noise (e.g. sensor noise). The sensor noise floor isassumed to be a mix of many independent causes and hence normal by thecentral limit theorem. In the phasor domain, this manifests as acircularly-symmetric complex normal distribution as given above.Unfortunately sensor noise floor is typically specified in dB(A), whichmeans only an upper bound is known for each frequency. In our model thenoise is assumed to be generated after conversion from acoustical to theelectrical domain.Step b

Convert the model from the previous step to a multivariate complexlog-normal distribution:

Gain Approximation:μ_(m)

√{square root over (1+σ²)}  (6)

$\begin{matrix}{\mu_{g}\overset{def}{=}{\frac{\sigma^{2}}{2\;\mu_{m}^{2}} - {\log\;\mu_{m}} + {\log{{s \cdot {w \odot x}}}}}} & (7) \\{\sigma_{g}\overset{def}{=}\;\frac{\sigma}{\mu_{m}}} & (8)\end{matrix}$log|o||s, w, x˜

(log|o|; μ _(g), σ_(g) ² I)  (9)

For phase:σ_(p)

σ  (10)∠o|s, w, x˜

(∠o; ∠[s·w⊙x], σ_(p) ² I)  (11)Combining Gain and Phase Gives:o|s, w, x˜

(o; μ _(g) +i∠[s·w└x], [σ_(g) ² +iσ _(p) ²]I)  (12)The approximations of step b are accurate when the signal to noise ratioof the observation is 7 dB or better.Step c

Model the uncertainties in the source signal and waveguide transferfunction as another multivariate complex log-normal distribution.

Source Signal:s˜

(s; μ _(s), Σ_(s))  (13)Waveguide Transfer Function:w˜

(w; μ _(w); Σ_(w))  (14)Combined Distribution:μ_(t)

μ_(s)+μ_(w)  (15)Σ_(t)

Σ_(s)+Σ_(w)  (16)t

s·w˜

(t; μ _(t), Σ_(t))  (17)A suitable waveguide/transmission medium has strong correlations betweenthe paths from source to each of the sensors in the waveguide. Thisamounts to the requirement that any undesired deviations in the behaviorof the waveguide, for example as caused by environmental factors such astemperature or the operator who is performing the measurement, shouldapply equally to each path in the transfer function w. When a suitablewaveguide is used, assumptions about the test source can be (very) mild.In the simplest case, the waveguide is the free field, e.g. whenpositioning the source and microphone array at known positions in ananechoic room. When such a controlled environment is not available, itcan be advantageous to couple the source more tightly to the array usinga waveguide.

In practice this distribution often has location parameter set to zero,and scale based on the assumptions about source and waveguide. The scaleparameters also encode the strong correlations that are typicallypresent in a suitable waveguide. Unknown properties (e.g. the phase ofthe source) can be modelled as infinite (or in a practicalimplementation: very large) scale parameters. It is assumed that thegain and phase uncertainties here are independent. For example, a lackof information can be incorporated in the model as follows. For thephase, set the mean to 0 and the variance to π or more. This causes theprior on the source phase to become approximately uniform due to thetails of the normal distribution of phase wrapping around in thecorresponding exponential. For the gain, the mean and variance can beset so that the prior covers the full dynamic range of the source.

Step d

Multiply the distributions of (b) and (c), then integrate out thenuisance variable t. This produces the multivariate complex log-normallikelihood distribution:p(o, t|x)=p(o|t, x)p(t)  (18)p(o|x)=∫p(o, t|x)dt  (19)It is easiest to do this for the gain and phase parts separately. Forthe gain:

$\begin{matrix}\begin{matrix}{{p\left( {{\log{o}},{{\log{t}}❘x}} \right)} = {{p\left( {{{\log{o}}❘{\log{t}}},x} \right)} \cdot {p\left( {\log{t}} \right)}}} \\{= {{\mathcal{N}\left( {{{\log{o}};\mu_{g}},{\sigma_{g}^{2}I}} \right)} \cdot {\mathcal{N}\left( {{{\log{t}};{\Re\;\mu_{t}}},{\Re\sum\limits_{t}}} \right)}}} \\{= {{\mathcal{N}\left( {{{\log{o}};{\frac{\sigma^{2}}{2\mu_{m}^{2}} - {\log\;\mu_{m}} + {\log{x}} + {\Re\mu}_{t}}},{{\sigma_{g}^{2}I} + {\Re\sum\limits_{t}}}} \right)} \cdot}} \\{{\mathcal{N}\left( {{\log{t}};\ldots}\mspace{11mu} \right)},}\end{matrix} & (20)\end{matrix}$Where the last equality follows by considering both factors as functionsof log|t| and applying known results. Because log|t| does not appear inthe first factor, integrating it out is trivial:

$\begin{matrix}{{p\left( {{\log{o}}❘x} \right)} = {{\int{{p\left( {{\log{o}},{{\log{t}}❘x}} \right)}{dt}}} = {{\mathcal{N}\left( {{{\log{o}};{\frac{\sigma^{2}}{2\;\mu_{m}^{2}} - {\log\;\mu_{m}} + {\log{x}} + {\Re\mu}_{t}}},{{\sigma_{g}^{2}I} + {\Re\sum\limits_{t}}}} \right)}.}}} & (21)\end{matrix}$Similarly, for the phase:p(∠o|x)=

(∠o; ∠x+

μ _(t), σ_(p) ² I+

Σ _(t))  (22)Combining:

$\begin{matrix}{\mu_{l}\overset{def}{=}{\frac{\sigma^{2}}{2\mu_{m}^{2}} - {\log\;\mu_{m}} + {\log{x}} + {i\;\angle\; x} + \mu_{t}}} & (23) \\{\sum\limits_{l}{\overset{def}{=}{{\left( {\sigma_{g}^{2} + {i\;\sigma_{p}^{2}}} \right)I} + \sum\limits_{t}}}} & (24) \\{o❘{\left. x \right.\sim{{\mathcal{C}\mathcal{L}\mathcal{N}}\left( {{o;\mu_{l}},\sum\limits_{l}} \right)}}} & (25)\end{matrix}$Step e

Model the prior sensor calibration state, again as a multivariatecomplex log-normal distribution:x˜

(x; μ _(x), Σ_(x))  (26)The location parameters of this distribution are set to the nominalsensor sensitivity and phase offset. The scale parameters are based onthe tolerances in sensor sensitivity and phase offset. Again, unknowntolerances can be modelled as infinite (or in a practicalimplementation: very large) scale parameters. It is assumed that thegain and phase tolerances of the sensors are independent.Step f

Bayes' rule is applied to get the posterior:x|o˜

(o; μ ₁, Σ₁)

(x; μ _(x)Σ_(x))/Z,  (27)where Z is a normalization constant. The likelihood (d) and prior (e)distributions can now be separated into real and imaginary parts (due tothe circularly-symmetric property (under (a)) and the independenceproperty (under (c) and (e)). These parts correspond to the gain andphase of the sensor, respectively. The gain and phase parts areprocessed separately.Gain

$\begin{matrix}\begin{matrix}{{p\left( {{\log{x}}❘{\log{o}}} \right)} = {{\mathcal{N}\left( {{{\log{o}};{\Re\mu}_{l}},{\Re\sum\limits_{l}}} \right)}{{\mathcal{N}\left( {{{\log{x}}❘{\Re\mu}_{x}},{\Re\sum\limits_{x}}} \right)}/Z_{1}}}} \\{= {\mathcal{N}\left( {{{\log{x}};{{- \frac{\sigma^{2}}{2\mu_{m}^{2}}} + {\log\;\mu_{m}} + {\log{o}} - {\Re\mu}_{t}}},{\Re\sum\limits_{l}}} \right)}} \\{{\mathcal{N}\left( {{{\log{x}}❘{\Re\mu}_{x}},{\Re\sum\limits_{x}}} \right)}/Z_{1}} \\{{= {\mathcal{N}\left( {{{\log{x}};{\Re\mu}_{c}},{\Re\sum\limits_{c}}} \right)}},}\end{matrix} & (28)\end{matrix}$where

μ_(c) and

Σ_(c) follow from the Wiener filter:

$\begin{matrix}{\mspace{85mu}{{{\Re\sum\limits_{c}}\overset{def}{=}\left\lbrack {\left( {\Re\sum\limits_{l}} \right)^{- 1} + \left( {\Re\sum\limits_{x}} \right)^{- 1}} \right\rbrack^{- 1}},}} & (29) \\{{\Re\mu}_{c}\overset{def}{=}{\left( {\Re\sum\limits_{c}} \right){\quad{\left\lbrack {{\left( {\Re\sum\limits_{l}} \right)^{- 1}\left( {{- \frac{\sigma^{2}}{2\mu_{m}^{2}}} + {\log\;\mu_{m}} + {\log{o}} - {\Re\mu}_{t}} \right)} + {\left( {\Re\sum\limits_{x}} \right)^{- 1}\left( {\Re\mu}_{x} \right)}} \right\rbrack.}}}} & (30)\end{matrix}$Phase

Similarly, for the phase:

$\begin{matrix}\begin{matrix}{{p\left( {{\angle\; x}❘{\angle\; o}} \right)} = {{\mathcal{N}\left( {{{\angle\; o};{\mathcal{J}\mu}_{l}},{\mathcal{J}\sum\limits_{l}}} \right)}{{\mathcal{N}\left( {{{\angle\; x}❘{\mathcal{J}\mu}_{x}},{\mathcal{J}\sum\limits_{x}}} \right)}/Z_{2}}}} \\{= {{\mathcal{N}\left( {{{\angle\; x};{{\angle\; o} - {\mathcal{J}\mu}_{t}}},{\mathcal{J}\sum\limits_{l}}} \right)}{{\mathcal{N}\left( {{{\angle\; x}❘{\mathcal{J}\mu}_{x}},{\mathcal{J}\sum\limits_{x}}} \right)}/Z_{2}}}} \\{{= {{\gamma\left( {\angle\; o} \right)}{{\mathcal{N}\left( {{{\angle\; x};{\mathcal{J}\mu}_{c}},{\mathcal{J}\sum\limits_{c}}} \right)}/Z_{2}}}},}\end{matrix} & (31)\end{matrix}$where

μ_(c) and

Σ_(c) again follow from the Wiener filter:

Σ _(c)

[(

Σ₁)⁻¹+(

Σ_(x))⁻¹]⁻¹,  (32)

μ _(c)=(

Σ _(c))[(

Σ₁)⁻¹(∠o−

μ _(t))+(

Σ_(x))⁻¹(

μ_(x))],  (33)γ(∠o)

(∠o;

μ _(t)+

μ_(x),

Σ₁+

Σ_(x)).  (34)The normalization constant γ/Z₂ has been retained explicitly, because acomplication arises: the (unwrapped) angle ∠o is unobservable and musthence be replaced with the (measured) angle

o+2πk, where k denotes the phase ambiguity. Hence:

$\begin{matrix}{{p\left( {{\angle\; x}❘{\measuredangle\; o}} \right)} = {\frac{1}{z_{3}}{\sum\limits_{k \in {\mathbb{Z}}^{N}}{{\gamma\left( {{\measuredangle\; o} + {2\pi\; k}} \right)}{{\mathcal{N}\left( {{{\angle\; x};{\mathcal{J}\;{{\overset{\_}{\mu}}_{c}\left( {{\measuredangle\; o} + {2\;\pi\; k}} \right)}}},{\mathcal{J}\;{\sum\limits^{\_}}_{c}}} \right)}.}}}}} & (35)\end{matrix}$Note that

μ _(c) now depends on the phase ambiguity k.

This expression can be interpreted as an infinite weighted mixture ofWiener filters (or an infinite weighted sum of normal distributions).For any practical implementation, it is required to truncate thisinfinite mixture and keep only the terms with dominant weights.

Unfortunately, selecting these terms amounts to the NP-hard ClosestLattice Vector problem. Moreover, existing heuristic methods, e.g. asdescribed by Morelande in (IEEE Conference on Acoustics, Speech andSignal Processing, 2008, pp3441-3444), do not cope well with the strongcorrelations encoded in step (c). Therefore, we introduce a novelselection algorithm called shotgun unwrapping. The result of applyingshotgun unwrapping (described in detail below) to truncate the summationand keep only the terms for which γ is significant is a pruned set{circumflex over (k)}∈

.

Finally, the posterior calibration mean and covariance are extracted bysummarizing the mixture distribution (using expected value identities):

μ_(c)

γ(

o+2π{circumflex over (k)})

μ _(c)(

o+2π{circumflex over (k)})/

γ(

o+2π{circumflex over (k)})  (36)

Σ_(c)

Σ _(c)+

γ(

o+2π{circumflex over (k)})[

μ _(c)(

o+2π{circumflex over (k)})]²/

γ(

o+2π{circumflex over (k)})−

μ_(c) ²  (37)Alternatively, the posterior mean and covariance can be evaluated forthe circular variableexp i∠x|o.  (38)This can give better results in practice, but has more complicated‘circular moment’ expressions. For the circular mean:

μ_(c,circular)

∠{

γ(

o+2π{circumflex over (k)})exp[i

μ _(c)(

o+2π{circumflex over (k)})]}  (38a)A full circular covariance matrix could be constructed using circularcorrelation coefficients, but in practice the circular standarddeviations of the individual microphones are sufficient to provide errorbars on the calibration:

$\begin{matrix}{\sigma_{c,{circular}}\overset{def}{=}\sqrt{{- 2}\;\ln{{\sum\limits_{\hat{k} \in \mathcal{K}}{{\gamma\left( {{\measuredangle\; o} + {2\;\pi\;\hat{k}}} \right)}{{\exp\left\lbrack {{i\;\mathcal{J}\;{{\overset{\_}{\mu}}_{c}\left( {{\measuredangle\; o} + {2\;\pi\;\hat{k}}} \right)}} - {\frac{1}{2}{diag}\;\mathcal{J}\sum\limits_{c}^{\_}}} \right\rbrack}/{\sum\limits_{\hat{k} = \mathcal{K}}{\gamma\left( {{\measuredangle\; o} + {2\;\pi\;\hat{k}}} \right)}}}}}}}} & \left( {38b} \right)\end{matrix}$Shotgun Unwrapping

An expression for the weights as function of the summation index can bederived as follows:

$\begin{matrix}{{\gamma(k)} = {\mathcal{N}\left( {{{{\measuredangle\; o} + {2\;\pi\; k}};{{\mathcal{J}\mu}_{t} + {\mathcal{J}\mu}_{x}}},{\mathcal{J}{\sum\limits_{l}{{+ \mathcal{J}}\sum\limits_{x}}}}} \right)}} \\{= {\mathcal{N}\left( {{{2\;\pi\; k};{{\mathcal{J}\;\mu_{t}} + {\mathcal{J}\;\mu_{x}} - {\measuredangle\; o}}},{\mathcal{J}{\sum\limits_{l}{{+ \mathcal{J}}\sum\limits_{x}}}}} \right)}} \\{= {\mathcal{N}\left( {{k;\frac{{\mathcal{J}\mu}_{t} + {\mathcal{J}\mu}_{x} - {\measuredangle\; o}}{2\;\pi}},\frac{\mathcal{J}{\sum\limits_{l}{{+ \mathcal{J}}\sum\limits_{x}}}}{4\;\pi^{2}}} \right)}} \\{\overset{def}{=}{\mathcal{N}\left( {{k;\mu_{u}},\sum\limits_{u}} \right)}}\end{matrix}$A large amount of points (‘pellets’) is drawn from the (continuous)random variable:k ₁˜

(k ₁; μ_(u), Σ_(u)),  (40)which are denoted by the set

₁ of size K₁. This can be done efficiently by first sampling from astandard multivariate normal distribution and then coloring and addingthe mean.

The pellets are rounded to their nearest integer value and duplicatesare removed. Each pellet now corresponds to a discrete value of k.Denote this as the set

₂ of size K₂.

Finally, the pellets are pruned by removing those with small weights.This can be done by sorting the pellets by Mahalanobis distance for eachpellet, defined as:M(k)

∥k−μ _(u)∥_(Σ) _(u) =√{square root over ((k−μ _(u))^(H)Σ_(u) ⁻¹(k−μ_(u)))},  (41)discarding all pellets with distance larger than a threshold (e.g. 0.95equiprobability curve of k₁), and finally returning the K shortest ones,where K is a practical upper limit for the amount of terms to beconsidered. The K selected pellets are denoted as {circumflex over (k)}∈

.

With high probability, the algorithm returns a good set of phaseunwrappings; the likelihood that a pellet ends up in the EuclideanVoronoi region of a lattice point is similar to the likelihood of thelattice point itself.

TABLE 1 Notation. Notation Meaning a~f(a) a has probability densityfunction f(a) a~f(a) a approximately has probability density functionf(a) p(a) Probability density function of a ⊙ Elementwise (Hadamard)product |a| Elementwise absolute value <a Elementwise argument a|bRandom variable a, given a realization of b

TABLE 2 Input parameters Known (input) parameters Symbol Domain MeaningUnit N

 ₊ Number of sensors — σ

Sensor noise floor Pa μ_(s)

Expected signal from the source — Σ_(s)

Uncertainty (variance) in signal from the — source μ_(w)

 ^(N) Nominal waveguide transfer functions — Σ_(w)

 ^(N×N) Uncertainty and correlations in waveguide — transfer functionsμ_(x)

 ^(N) Nominal sensor gain ( 

 ) and phase ( 

 ) — before calibration Σ_(x)

 ^(N×N) Tolerances in sensor gain ( 

 ) and phase — ( 

 ) before calibration

TABLE 3 Output parameters Unknown (output) parameters Symbol DomainMeaning Unit μ_(c)

 ^(N) Nominal sensor gain ( 

 ) and phase ( 

 ) — after calibration Σ_(c)

 ^(N×N) Tolerances in sensor gain ( 

 ) and phase — ( 

 ) after calibration

TABLE 4 Intermediate variables for the regression Intermediate variables(regression) Symbol Domain Meaning Unit μ_(m)

Intermediate mean in approximation of Pa observation gain μ_(g)

 ^(N) Mean of log-normal approximation of — observation gain σ_(g)

Std. dev. of log-normal approximation — of observation gain σ_(p)

Std. dev. of log-normal approximation — of observation phase μ_(t)

 ^(N) Mean of source-waveguide product — Σ_(t)

 ^(N×N) Covariance of source-waveguide product — Z, Z _(k)

Unimportant normalization constants —

 ū_(c) ({circumflex over (k)})

 ^(N) Nominal phase after calibration for — specific phase unwrapping

 Σ _(c)

 ^(N×N) Tolerances in phase after calibration for — specific unwrapping

TABLE 5 Intermediate variables for the phase unwrapping Intermediatevariables (phase unwrapping) Symbol Domain Meaning Unit k

 ^(N) Phase ambiguity — {circumflex over (k)}

Resolved phase ambiguity — μ_(u)

 ^(N) Mean for shotgun unwrapping — Σ_(u)

 ^(N×N) Covariance for shotgun unwrapping — k₁

 ^(N) Continuous phase ambiguity for — shotgun unwrapping

 ₁ { 

 ^(N)}^(K) ₁ The set of K₁ ∈ Z₊ samples from k₁ —

 ₂ { 

 ^(N)}^(K) ₂ The set of K₂ ∈ Z₊ ≤ K₁ unique and — rounded samples from

 ₁

{ 

 ^(N)}^(K) The set of K ∈ Z₊ ≤ K₂ most — dominant samples from

 ₂

TABLE 6 Unobserved random variables Unobserved (latent) random variablesSymbol Domain Meaning Unit x

 ^(N) Calibration state — s

Source signal Pa w

 ^(N) Waveguide transfer functions — t

 ^(N) Product of source signal and Pa waveguide transfer functions —

TABLE 7 Observed random variables Observed (measured) random variablesSymbol Domain Meaning Unit o

 ^(N) Observation Pa

The invention claimed is:
 1. A method of calibrating gains and/or phasesof elements of an array of N acoustic microphones, the methodcomprising: providing an acoustic source; providing an estimate of atransfer function from the acoustic source to the elements of the arrayof N acoustic microphones; performing one or more measurements ofacoustic signals received at the elements of the array of N acousticmicrophones when the acoustic source is operating; performing Bayesianinference of gains and/or phases of the array of N acoustic microphonesbased at least on the one or more measurements and on the estimate ofthe transfer function.
 2. The method of claim 1, wherein a posteriorphase probability distribution of the Bayesian inference is an infiniteweighted sum of normal distributions, each normal distribution having acorresponding weight γ(k), where k is an N-dimensional vector ofintegers; wherein a phase unwrapping of the Bayesian inference isperformed by sampling a probability distribution of γ(k) to provide ak-set and selecting the K best values from the k-set, wherein K is apredetermined integer.
 3. The method of claim 2, wherein sampling aprobability distribution of the weights γ(k) to provide a k-set andselecting the K best values from the k-set comprises: sampling from acontinuous probability distribution of γ(k) to provide an initial k-set

₁; rounding elements of the initial k-set

₁to the nearest integers and eliminating any resulting duplicates toprovide a discretized k-set

₂; evaluating distances of each element of

₂ from a mean of the probability distribution of γ(k); selecting the Kelements of

₂ having the shortest distances as the K best values.
 4. The method ofclaim 3 wherein the selecting the K elements of

₂ having the shortest distances comprises removing elements of

₂ having distances greater than a predetermined threshold prior toselecting the K best weights.
 5. The method of claim 3, wherein theprobability distribution of γ(k) has a mean μ_(u) and a covariancematrix Σ_(u), and wherein the evaluating distances M(k) comprisescalculating M(k) =√{square root over ((k−μ_(u))^(H)Σ_(u) ⁻¹(k−μ_(u)))} .6. The method of claim 1, wherein amplitude and/or phase of the acousticsource are assumed to be drawn from a predetermined source probabilitydistribution.
 7. The method of claim 1, wherein the transfer functionsare determined by an acoustic waveguide network configured to couple theacoustic source to the array of acoustic microphones.
 8. The method ofclaim 7, wherein the acoustic waveguide network includes a source portcorresponding to the acoustic source, and array ports, each array portcorresponding to a corresponding one of the elements of the array ofacoustic microphones.
 9. The method of claim 1, wherein the acousticsource is an uncalibrated acoustic source.
 10. The method of claim 9,wherein the acoustic source is part of a mobile electronic device. 11.The method of claim 1, wherein the acoustic source comprises an acousticcalibrator or pistonphone.
 12. The method of claim 1, further comprisingusing an auxiliary reference microphone to provide a traceablecalibration of the array of N acoustic microphones.
 13. The method ofclaim 1, wherein the Bayesian inference is further based on informativeprior estimates of gains and/or phases of the array of N acousticmicrophones.
 14. The method of claim 13, wherein the informative priorestimates of gains and/or phases of the array of N acoustic microphonesare derived from manufacturer specifications of the array of N acousticmicrophones.